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ABSTRACT 

Measurements of the cosmic ray p/p ratio are compared to predictions from an 
inhomogeneous disk-diffusion model of p production and propagation within the 
Galaxy, combined with a calculation of the modulation of the interstellar cosmic 
ray spectra as the particles propagate through the heliosphere to the Earth. The 
predictions agree with the observed p/p spectrum. Adding a finite p lifetime to 
the model, we obtain the limit Tp > 0.8 Myr (90% C.L.). 
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1. Introduction 

In recent years the presence of antiprotons (p's) in the cosmic ray (CR) flux incident 
upon the Earth has been firmly established by a series of balloon-borne experiments (Golden 
et al. 1979; Bogomolov et al. 1979; Bogomolov et al. 1987; Bogomolov et al. 1990; Hof 
et al. 1996; Mitchell et al. 1996; Moiseev et al. 1997; Yoshimura et al. 1995; Boezio et 
al. 1997; Matsunaga et al. 1998). The measurements are summarized in Table 1. The 
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observed CR p/p ratio has been shown to be in approximate agreement with predictions 
based on the Leaky Box Model (LBM) (Stephens 1981; Stephens k Golden 1987; Webber & 
Potgieter 1989; Gaisser & Schaefer 1992), which assumes that the p's originate from proton 
interactions in the interstellar medium (ISM). The p's then propagate within the Galaxy 
until they "leak out" by diffusion and convection with the characteristic CR Galactic escape 
time T ~ 10 million years (Myr) (Webber & Potgieter 1989; Gaisser & Schaefer 1992). If 
the p lifetime Tp is not long compared to T the predicted p/p spectrum will be modified. 
The agreement of the LBM predictions with the observed p/p spectrum has therefore been 
used to argue that Tp > 10 Myr (Golden et al. 1979; Bogomolov et al. 1979; Stephens 1981; 
Stephens & Golden 1987). This estimated limit is based on early CR p data and does not 
take into account the reduction of the p decay rate due to time dilation, the effect of the 
heliosphere on the observed p/p spectrum, or the systematic uncertainties associated with the 
predictions. In this paper we compare recent CR data with the predictions of an improved 
LBM extended to permit a finite Tp. Heliospheric corrections and systematic uncertainties 
are taken into account. Assuming a stable p, we find excellent agreement between our 
predictions and the CR observations. Allowing the p to decay, we obtain a lower limit on Tp 
which is significantly more stringent than current laboratory bounds obtained from searches 
for p decay in ion traps (Gabrielse et al. 1996) and storage rings (Geer et al. 1994; Hu et 
al. 1998). The analysis presented in this paper improves on our earlier analysis (Geer & 
Kennedy 1998) by including new data from Matsunaga et al. (1998). 

CPT invariance requires Tp — Tp., where the proton lifetime Tp is known to exceed 0(10^^) 
yr (Caso et al. 1998). Although there is no compeUing theoretical motivation to suspect a vi- 
olation of CPT invariance, and hence a short p lifetime, it should be noted that string theories 
can accommodate CPT violation. Consider a mass-dimension-n CPT-violating quantum 
field operator suppressed by the characteristic scale mx, with n > 4. Dimensional analysis 
provides the estimate rripTp ~ [mp/mx]^"~^, yielding mx/rrip ~ [4.5 x lO^^-rp/lO Myr]-*^/*^^""^) 
(Kennedy 1999). For a given lower limit on r^, the implied lower limit on nix is most strin- 
gent for n — b. Note that if mx is at the Planck scale (1.2 x 10^^ GeV/c^) and n — the 
expected Tp would be ~ 10 Myr. Hence, a search for p decay with a lifetime approaching 
10 Myr provides a test for CPT violation well beyond the scale accessible at high energy 
colliders. Finally, since the antiproton is the only long lived antiparticle that could decay 
into other known particles without violating charge conservation, a search for a modification 
of the CR p spectrum by p decay provides a unique test of the stability of antimatter. 
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2. Galactic Production and Propagation 

In the LBM the ISM p's are produced by the interactions of CR p's (Stephens 1981; 
Stephens & Golden 1987; Webber & Potgieter 1989; Gaisser & Schaefer 1992): pNz pX, 
where Nz is a nucleus of charge Z, and X is anything. Our calculations use the elemental 
ISM abundances given in Webber & Potgieter (1989) and Gaisser & Schaefer (1992), and 
the measured cross sections ior Z — 1 (the dominant contribution) given in Stephens (1981). 
For Z > 1, we have used the "wounded nucleon" results also used by Gaisser & Schaefer 
(1992). 

The p's are assumed to propagate within the Galaxy until they are lost either by leakage 
into intergalactic space or by pp annihilation. The dominant loss process is leakage. Energy 
loss is included in the model. Gomparisons between predicted and observed GR isotopic 
abundances imply that reacceleration can be neglected (Webber et al. 1992). Reacceleration 
is not included in our calculation. The LBM incorporates space-dependent diffusive and 
convective CR leakage into a single escape term with a single characteristic escape time. 
This picture is known to be too simple. A more accurate picture of Galactic GRs (Webber 
et al. 1992) is provided by a two-zone diffusion-convection model which we refer to as the 
inhomogeneous leaky disk model (ILDM). The model is based on a /i ~ 100 pc-thick matter 
disk (radius i? ~ 20 kpc), with a halo extending above and below the disk by L ~ ±3 kpc. 
The GRs diffuse off the disk into the halo by magnetic turbulence. It has been shown that 
convection driven by the Galactic wind can be neglected (Bloemen 1993). The ILDM disk 
storage time ~ where the diffusion transverse to the Galactic disk /C — (6 ±4) x 10^'^ 

m^ sec~^ (Bloemen 1989; Webber et al. 1992). The ILDM with these parameters (h, R, L, K,) 
has been numerically solved by Ghardonnet et al. (1996) for the secondary GR p spectrum. 
Although the ILDM provides a more realistic physical description of GR propagation than 
the LBM, Ghardonnet et al. found that the fluxes in the Galactic disk predicted by the ILDM 
are reproduced within uncertainties by a modified LBM (MLBM) in which the momentum- 
dependent escape time T{P) is given by: T{P) = (13 Myr) [l+cP/(3 GeV)]-°-^. This escape 
time is longer than that used by Webber & Potgieter (1989) and Gaisser & Schaefer (1992) 
and, since the p flux is proportional to T{P), raises the predicted p flux and continues the 
trend of the last two decades toward longer GR escape times/lengths, as discussed by these 
authors. 

The uncertainties on the parameters of the MLBM result in uncertainties on the normal- 
ization of the predicted p/p ratio but, to a good approximation, do not introduce significant 
uncertainties in the shape of the predicted spectrum. Uncertainties on four parameters must 
be considered: (i) the storage time (±67% ) (Webber et al. 1992), (ii) the ISM primary p flux 
(±35% ) (Gaisser & Schaefer 1992), (iii) the p production cross section (±10% ) (Stephens 
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1981; Gaisser & Schaefer 1992), and (iv) the composition of the ISM, which introduces an 
uncertainty of < 6% on the predicted p flux (Gaisser & Schaefer 1992). We neglect the last 
of these uncertainties since it is relatively small. Within the quoted fractional uncertain- 
ties on the other three parameters wc treat all values as being a priori equally likely. Note 
that the predicted p/p ratio is approximately proportional to each of the parameters under 
consideration. 



3. Heliospheric Modulation 

The solid curves in Figure 1 show the MLBM p/p spectra for the parameter choices 
that result in the largest and smallest p/p predictions. The predicted ISM spectrum does 
not give a good description of the observed distribution at the top of the atmosphere. Good 
agreement is not expected because the CR spectra observed at the Earth are modulated as 
the particles propagate through the heliosphere, which consists of the solar magnetic field 
B and solar wind (Stix 1989; Encrenaz et al. 1990; Longair 1992). The wind, assumed 
to blow radially outwards, has equatorial speed Vw — 400 km sec~^. Ulysses has found a 
latitude-dependent Vw{9). At high latitudes (|^-90°| > 20°), Vw ^ 750 km sec"^ (Smith et 
al. 1995; Marsden et al. 1996; Ulysses 1999). The latitude-dependent solar rotation flQ{6) 
(Howard 1984) twists the field lines into a Parker spiral. The smoothed heliomagnetic field 
(i?0 ~ 5 nT at the Earth's orbit) declines as it changes from radial at the Sun to azimuthal 
in the outer solar system. The heliomagnetic polarity sgn{A) is opposite in northern and 
southern solar hemispheres and switches sign somewhat after sunspot maximum. Note that 
^4 > in the 1990-99 epoch. The regions of opposite magnetic polarity are separated 
by an approximately equatorial, unstable neutral current sheet. The sheet is spiraled and 
wavy; its waviness is measured by the tilt angle a, which relaxes from ~ 50° at polarity 
reversal to ~ 10° just before reversal. The wavy sheet locus in heliographic co-ordinates is 
cos^ + sina • sm[(p + flQr/Vw] — (Jokipii & Kopriva 1979; Encrenaz et al. 1990; NSO 
1999). 

Cosmic rays enter the hehosphere on baUistic trajectories and are then transported by 
drift and diffusion (Jokipii & Kopriva 1979; Riskcn 1984). The transport model includes 
energy loss as CRs perform work against the wind, CR solar wind convection, drift of CRs 
across inhomogeneous field lines, and diffusive scattering by field turbulence. Particles with 
qA > {qA < 0) drift in along a polar (sheet) route. Since the polar Vyy is approximately 
twice the equatorial Vw, particles with qA > find it more difficult to drift in along a polar 
route than do qA < particles along a sheet route. In the following, we restrict our analysis 
to CR p's and p's with kinetic energies > 500 MeV. At these energies, (a) the wind and 
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magnetic drift (Isenberg & Jokipii 1979, 1981) dominate the transport; and (b) the diffusion 
is charge-symmetric and, for a given rigidity, modulates the p's and p's by the same factor, 
and hence cancels in the ratio of fluxes (see Ulysses e and p scans: Smith et al. 1995; 
Marsden et al. 1996; Fisk et al. 1998). 

We compute the modulation of the CR fluxes by the method of characteristics, numer- 
ically implemented in three dimensions by a combination of Rungc-Kutta and Richardson- 
Burlich-Stoer techniques (Press et al. 1992). The calculation is based on the hcliospheric 
transport model of Jokipii & Kopriva (1979), Jokipii & Thomas (1981), and Jokipii & Davila 
(1981), but that model has been updated with Pioneer, Voyager, Helios, IMP and Ulysses 
hcliospheric measurements (Smith et al. 1995; Marsden et al. 1996; NSSDC 1999; SEC 
1999). The improvements include the latitude-dependent Vw{0) and ^q{9). The calculation 
is simplified by ignoring diffusion where its effects on the p/p ratio are small, which for 
our energy range means everywhere except across the sheet (Jokipii & Thomas 1981). The 
effects of diffusion away from the sheet are essential for particles with energies below the 
range we are considering. For particles with speed v the across-sheet diffusion coefficient k± 
= [(2 - 3) X 10^^ m^/sec][B(s/B{r)]{P/GeVy/^{v/c) is used (Smith et al. 1995; Marsden et 
al. 1996). The implied field fiuctuation spectral index is the Kolmogorov value of 5/3. The 
modulated fluxes at Earth depend only weakly on k± in our energy range, for this and larger 
values of k±. 

4. Results 

We use the data sets recorded by the MASS91, IMAX, BESS, and CAPRICE exper- 
iments (Table 1). These data were recorded in the period 1991-1995, corresponding to a 
well-behaved part of the solar cycle for which the hcliospheric modulation corrections can 
be confidently calculated. To explore the dependence of the predicted spectrum on the hc- 
liospheric parameters (equatorial Vw, polar Vw, and B^) we have computed the modulated 
spectra for 11 parameter sets (Fl - Fll) that span the range of acceptable parameter values 
(Table 2). Using the central parameter values for our MLBM, and assuming a stable an- 
tiproton, the predicted modulated p/p spectra for a fixed time in the solar cycle (July 1995) 
are shown in Figure 2 for each of the 11 hcliospheric parameter sets. Figure 2 also shows the 
epoch-corrected measured CR spectra, obtained by multiplying each measurement by the 
factor / = i?(July 1995)/i?(t), where R{t) is the predicted p/p ratio at time t. The factors 
/, which are shown in Table 1 and have been computed using the F6 parameters, vary by up 
to ±0.06 with the parameter set choice. The MLBM predicted p/p spectra give an excellent 
description of the measurements. There is no evidence for an unstable p. 
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To obtain a limit on Tp we add to the MLBM one additional loss mechanism, p decay. 
The results from maximum likelihood fits to the measurements are shown in Figure 3 as 
a function of the assumed Tp for the 11 heliospheric parameter sets. The fits, which take 
account of the Poisson statistical fluctuations on the number of observed events and the 
background subtraction for each data set (Table 1), also allow the normalization of the 
MLBM predictions to vary within the acceptable range (Figure 1). For a stable p, at the 95% 
C.L. all of the heliospheric parameter sets yield predictions that give reasonable descriptions 
of the observed p/p spectrum. Allowing for a finite Tp, the heliospheric parameter sets with 
larger wind speeds permit lower p lifetimes. This can be understood by noting that, as the 
wind speed increases, the predicted flux of polar-routed particles (protons in the present 
solar cycle) is depleted at kinetic energies < 10 GeV, which increases the predicted p/p 
ratio in this energy range. Antiproton decay would compensate for this distortion in the 
predicted spectrum. Hence our flts using the extreme parameter set Fll determine the 
limits on r-p. Under the assumption that there are no significant non-standard sources of 
cosmic ray antiprotons, we obtain the bounds: 

Tp > 0.8 Myr (90% C.L.) , 0.7 Myr (95% C.L.) , 0.5 Myr (99% C.L.) . (1) 

Our simple dimensional analysis suggests that if a dimension-n CPT-violating coupling 
results in antiproton decay, the mass scale mx at which this new physics takes place exceeds 
0(10^9) GeV/c2 [C(109) GeV/c2], forn = 5 (6). 

The limits (1) are significantly more stringent than those obtained from the most sensi- 
tive laboratory search for inclusive p decay (r^ > 3.4 months) (Gabrielse et al. 1996) or the 
most sensitive search for an exclusive p decay mode {Tp/B{p — > /i~7) > 0.05 Myr) (Hu et 
al. 1998). Note however that these limits are less restrictive than the estimate Tp > 10 Myr 
given in Caso et al. (1998), which takes no account of time dilation, systematic uncertainties, 
or propagation effects and which is therefore overconstraining. 

It is a pleasure to thank J. R. Jokipii (Lunar & Planetary Laboratory, Univ. Arizona) 
and E. J. Smith (Jet Propulsion Laboratory/CalTech) for their insights. This work was 
supported at Fermilab under grants U.S. DOE DE-AC02-76CH03000 and NASA NAG5- 
2788 and at the Univ. Florida, Institute for Fundamental Theory, under grant U.S. DOE 
DE-FG05-86ER40272. 
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Table 1: Summary of Cosmic Ray Antiproton Results 



Experiment 


Field 


Flight 




KE Range 


Cand- 


Back- 


Observed 


Prcdict- 




Pol.*^ 


Date 




(GeV) 


idates 


ground 


p/p Ratio 


tion'^ 


Golden et al. 1979^ 


+ 


June ig7g 


- 


5.6 - 12.5 


46 


18.3 


(5.2 ± 1.5) X 10-^ 


- 


Bogomolov et al. 1979^ 


+ 


1972-1977 


- 


2.0 - 5.0 


2 


- 


(6 ± 4) X 10"^ 


- 


Bogomolov et al. 1987''' 




1984-1985 




0.2 - 2.0 


1 




(Qtl") X 10-5 
(2.4t?:^) X 10-^ 




Bogomolov et al. 1990-I- 




1986-1988 




2.0 - 5.0 


3 






MASS91'* 


+ 


Sep. 1991 


1.1 


3.70-19.08 


11 


3.3 


(1.24l°j8) X 10-^ 


1.3 X 10-^ 


IMAX^* 


+ 


July 1992 




0.25 - 1.0 


3 


0.3 


{s.utli) X 10-5 


1.5 X 10-5 


IMAX*^ 


+ 


July 1992 


0.96 


1.0 - 2.6 


8 


1.9 


(5.36^|5) X 10-5 


6.5 X 10-5 


IMAX*^ 


+ 


July 1992 


1.1 


2.6 - 3.2 


5 


1.2 


(1.94+}-?) X 10-1 


1.1 X 10-^ 


BESsgs^-f 


+ 


July 1993 




0.20 - 0.60 


7 


~ 1.4 


(5.21^:^) X 10-6 


8.9 X 10-6 


CAPRICE^ 


+ 


Aug. 1994 


0.94 


0.6 - 2.0 


4 


1.5 


(2.5^?:^) X 10-5 


3.5 X 10-5 


CAPRICE^ 


+ 


Aug. 1994 


1.0 


2.0 - 3.2 


5 


1.3 


(i.gli:^) X 10-^ 


1.1 X 10-^ 


BESSgs''^* 


+ 


July 1995 


1.0 


0.175 - 0.3 


3 


0.17 


(7.81^1) X 10-6 




BESSgs''^* 


+ 


July 1995 


1.0 


0.3 - 0.5 


7 


0.78 


(TAttD X 10-6 


1.1 X 10-5 


BESSgS*"* 


+ 


July 1995 


1.0 


0.5 - 0.7 


7 


1.4 


(7.7l|^) X 10-6 


5.5 X 10-6 


BESSgs''* 


+ 


July 1995 


1.0 


0.7 - 1.0 


11 


2.8 


(1.01+5-7) X 10-5 


1.3 X 10-5 


BESSgS*"* 


+ 


July 1995 


1.0 


1.0 - 1.4 


15 


3.5 




3.1 X 10-5 



^ Heliomagnetic field polarity sgn(74). 
^ Epoch correction factor; see text. 

MLBM prediction using F6 heliospheric parameters; see text. 
^ Hof et al. 1996. ^ Mitchell et al. 1996. ^ Moiseev et al. 1997. s Boezio et al. 1997. 
^ Matsunaga et al. 1998. 

Not shown in Fig. 1 or used in analysis. * Not used in analysis. * Statistical and systematic 
uncertainties on ratio added in quadrature. 
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Table 2: Hcliospheric Parameter Sets 



Set 


Equatorial Vw 
(km sec~^) 


Polar Vw 

(km sec""-*^) 


-6e 
(nT) 


Fl 


375 


700 


4.0 


F2 


380 


710 


4.1 


F3 


385 


720 


4.2 


F4 


390 


730 


4.3 


F5 


395 


740 


4.4 


F6 


400 


750 


4.5 


FT 


405 


760 


4.6 


F8 


410 


770 


4.7 


F9 


415 


780 


4.8 


FIO 


420 


790 


4.9 


Fll 


425 


800 


5.0 
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Kinetic Energy (GeV) 



Fig. 1. — Observed p/p ratio at the top of Earth atmosphere (see Table 1). The sohd curves 
show the upper and lower interstellar ratios predicted by the MLBM described in the text, 
without solar modulation. The broken curve shows the MLBM prediction with the same 
parameters used for the modulated predictions of Fig. 2. 
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Fig. 2. — Observed p/p spectrum (kinetic energy > 500 MeV) at the top of the atmosphere 
compared with the MLBM predictions (see broken curve on Fig. 1) after modulation using 
the hehospheric parameter sets indicated (see Table 2). The curves are predictions for the 
spectrum observed in July 1995. The data have been corrected to correspond to this epoch 
(see text). 
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Fig. 3. — Fit results as a function of the assumed Tp for the eleven heliospheric parameter 
sets (Fl - Fll) shown in Table 2. 



